#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>

#define MAX_POLY_DEGREE 100  // 最大多项式次数

// ========================= 多项式操作工具 =========================
// 多项式化简：将多项式矩阵（如[[1,-x],[1,-y]]）展开为单个多项式系数
// 输入：多项式因子数组（每个因子是[1, -x_i]，即 1 - x_i*x）
// 输出：化简后的多项式系数（从高次到低次），out_len返回系数长度
double* polynomial_simplify(double** factors, int factor_rows, int factor_cols, int* out_len) {
    // 初始多项式：1（系数[1]）
    double* res = (double*)malloc(sizeof(double));
    res[0] = 1.0;
    int res_len = 1;

    // 逐个因子相乘（(a0+a1x+...)(1 - xi x)）
    for (int i = 0; i < factor_rows; i++) {
        double* factor = factors[i];  // 因子：[1, -xi]
        int factor_len = factor_cols;

        // 新多项式长度 = 原长度 + 因子长度 - 1
        int new_len = res_len + factor_len - 1;
        double* new_res = (double*)calloc(new_len, sizeof(double));

        // 多项式乘法：res * factor
        for (int j = 0; j < res_len; j++) {
            for (int k = 0; k < factor_len; k++) {
                new_res[j + k] += res[j] * factor[k];
            }
        }

        // 更新结果
        free(res);
        res = new_res;
        res_len = new_len;
    }

    // 去除末尾零（可选，与Python逻辑对齐）
    while (res_len > 1 && fabs(res[res_len - 1]) < 1e-10) {
        res_len--;
    }

    *out_len = res_len;
    return res;
}

// 数组扩展（前面补0，对应Python的[0]*n + list）
double* array_extend_front(double* arr, int arr_len, int pad_len, int* out_len) {
    *out_len = pad_len + arr_len;
    double* new_arr = (double*)calloc(*out_len, sizeof(double));
    memcpy(new_arr + pad_len, arr, arr_len * sizeof(double));
    free(arr);
    return new_arr;
}

// ========================= 霍纳多项式结构体（模拟Python HornerPolynomial） =========================
typedef struct {
    double* coefficients;  // 多项式系数（从高次到低次，如x²+1 → [1,0,1]）
    int degree;            // 多项式次数（系数长度-1）
} HornerPolynomial;

// 创建霍纳多项式
HornerPolynomial* horner_poly_create(double* coeffs, int len) {
    HornerPolynomial* poly = (HornerPolynomial*)malloc(sizeof(HornerPolynomial));
    poly->degree = len - 1;
    poly->coefficients = (double*)malloc(len * sizeof(double));
    memcpy(poly->coefficients, coeffs, len * sizeof(double));
    return poly;
}

// 释放霍纳多项式内存
void horner_poly_destroy(HornerPolynomial* poly) {
    if (!poly) return;
    free(poly->coefficients);
    free(poly);
}

// 霍纳法则求值（核心功能）
double horner_poly_evaluate(const HornerPolynomial* poly, double x) {
    double result = 0.0;
    for (int i = 0; i <= poly->degree; i++) {
        result = result * x + poly->coefficients[i];
    }
    return result;
}

// 多项式转字符串（如[1,0,1] → "x**2+1"）
char* horner_poly_to_str(const HornerPolynomial* poly) {
    char* str = (char*)malloc(1024 * sizeof(char));
    str[0] = '\0';
    bool first = true;

    for (int i = 0; i <= poly->degree; i++) {
        double coeff = poly->coefficients[i];
        if (fabs(coeff) < 1e-10) continue;

        // 符号处理
        if (!first && coeff > 0) {
            strcat(str, "+");
        }
        first = false;

        // 系数部分
        int order = poly->degree - i;
        if (order == 0) {
            // 常数项
            sprintf(str + strlen(str), "%.1f", coeff);
        } else {
            if (fabs(coeff) != 1.0) {
                sprintf(str + strlen(str), "%.1f*", coeff);
            }
            // 次数部分
            strcat(str, "x");
            if (order > 1) {
                sprintf(str + strlen(str), "**%d", order);
            }
        }
    }

    // 空多项式返回"0"
    if (first) {
        strcpy(str, "0");
    }
    return str;
}

// ========================= 牛顿插值核心算法 =========================
// 获取单个基多项式（对应Python get_polynomial）
double* get_polynomial(double** data, int data_rows, int n, int* out_len) {
    int poly_len = data_rows;
    double* arr = (double*)calloc(poly_len, sizeof(double));

    if (n == 0) {
        arr[poly_len - 1] = 1.0;
        *out_len = poly_len;
        return arr;
    }

    // 构建多项式因子矩阵：[[1, -data[i][0]] for i in 0..n-1]
    double** factors = (double**)malloc(n * sizeof(double*));
    for (int i = 0; i < n; i++) {
        factors[i] = (double*)malloc(2 * sizeof(double));
        factors[i][0] = 1.0;
        factors[i][1] = -data[i][0];  // -x_i
    }

    // 多项式化简
    int simplify_len;
    double* p = polynomial_simplify(factors, n, 2, &simplify_len);

    // 释放因子矩阵
    for (int i = 0; i < n; i++) {
        free(factors[i]);
    }
    free(factors);

    // 前面补0，使长度等于data_rows
    int pad_len = data_rows - simplify_len;
    double* res = array_extend_front(p, simplify_len, pad_len, out_len);

    return res;
}

// 牛顿插值主函数（对应Python interpolate）
HornerPolynomial* newton_interpolate(double** data, int data_rows) {
    if (data_rows == 0) return NULL;

    // 步骤1：生成基多项式数组e
    double** e = (double**)malloc(data_rows * sizeof(double*));
    int* e_lens = (int*)malloc(data_rows * sizeof(int));  // 每个基多项式的长度
    for (int i = 0; i < data_rows; i++) {
        e_lens[i] = 0;
        e[i] = get_polynomial(data, data_rows, i, &e_lens[i]);
    }

    // 步骤2：计算差商表f
    double** f = (double**)malloc(data_rows * sizeof(double*));
    for (int i = 0; i < data_rows; i++) {
        f[i] = (double*)calloc(data_rows - i, sizeof(double));
    }

    // 初始化差商表第一列（f[i][0] = y_i）
    for (int i = 0; i < data_rows; i++) {
        f[i][0] = data[i][1];
    }

    // 计算高阶差商
    for (int j = 1; j < data_rows; j++) {
        for (int i = 0; i < data_rows - j; i++) {
            double xi = data[i + j][0];       // x_{i+j}
            double x0 = data[i][0];           // x_i
            f[i][j] = (f[i + 1][j - 1] - f[i][j - 1]) / (xi - x0);
        }
    }

    // 步骤3：求和得到最终多项式系数
    int poly_len = data_rows;
    double* polynomial = (double*)calloc(poly_len, sizeof(double));
    for (int i = 0; i < data_rows; i++) {
        double fi0 = f[0][i];  // 差商表第一行第i列
        double* ei = e[i];     // 第i个基多项式
        int ei_len = e_lens[i];

        // 同幂次系数相加
        for (int j = 0; j < poly_len; j++) {
            if (j < ei_len) {
                polynomial[j] += fi0 * ei[j];
            }
        }
    }

    // 步骤4：创建霍纳多项式
    HornerPolynomial* res_poly = horner_poly_create(polynomial, poly_len);

    // 释放临时内存
    free(polynomial);
    for (int i = 0; i < data_rows; i++) {
        free(e[i]);
        free(f[i]);
    }
    free(e);
    free(e_lens);
    free(f);

    return res_poly;
}

// ========================= 测试代码 =========================
int main() {
    // 测试数据：[[0,1],[2,5],[4,17]]
    int data_rows = 3;
    double** data = (double**)malloc(data_rows * sizeof(double*));
    data[0] = (double*)malloc(2 * sizeof(double));
    data[1] = (double*)malloc(2 * sizeof(double));
    data[2] = (double*)malloc(2 * sizeof(double));
    data[0][0] = 0.0; data[0][1] = 1.0;
    data[1][0] = 2.0; data[1][1] = 5.0;
    data[2][0] = 4.0; data[2][1] = 17.0;

    // 牛顿插值计算
    HornerPolynomial* poly = newton_interpolate(data, data_rows);
    if (!poly) {
        fprintf(stderr, "插值失败\n");
        return -1;
    }

    // 测试求值
    printf("poly(0) = %.1f, poly(2) = %.1f, poly(4) = %.1f\n",
           horner_poly_evaluate(poly, 0.0),
           horner_poly_evaluate(poly, 2.0),
           horner_poly_evaluate(poly, 4.0));

    printf("poly(1) = %.1f, poly(3) = %.1f, poly(5) = %.1f\n",
           horner_poly_evaluate(poly, 1.0),
           horner_poly_evaluate(poly, 3.0),
           horner_poly_evaluate(poly, 5.0));

    // 打印多项式字符串和系数
    char* poly_str = horner_poly_to_str(poly);
    printf("多项式表达式：%s\n", poly_str);
    printf("多项式系数：[");
    for (int i = 0; i <= poly->degree; i++) {
        printf("%.1f", poly->coefficients[i]);
        if (i != poly->degree) {
            printf(", ");
        }
    }
    printf("]\n");

    // 释放内存
    free(poly_str);
    horner_poly_destroy(poly);
    for (int i = 0; i < data_rows; i++) {
        free(data[i]);
    }
    free(data);

    return 0;
}